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Abstract — Interactions  of  the  detonation-product  gas,  shell-casing  fragments,  soil  ejecta  and  various 
other  debris  with  bulk  water  barriers  surrounding  the  explosive  have  been  demonstrated  to  have  a 
potentially  major  beneficial  effect  in  mitigation  cf  the  effects  of  an  explosion.  In  the  present  work, 
various  computational  methods  ranging  from  those  based  on  thermo-chemistry  of  the 
detonation/combustion  chemical  reactions  to  those  involving  transient,  nonlinear-dynamics  based 
mechanical  interactions  between  detonation  products,  air  and  water  are  used  to  better  understand  and 
quantify  the  beneficial  effects  of  various  potential  explosion-mitigation  mechanisms.  In  particular,  the 
absorption  of  the  detonation  energy  by  water,  water-aerosolization  induced  reduction  in  the  shock 
speed,  transfer  of  momentum  from  the  explosion  products  to  water  and  deceleration/suppression  of  the 
combustion  reactions  are  examined  computationally.  The  results  obtained  show  that  water  evaporation 
which  consumes  a  substantial  portion  of  the  detonation  energy  plays  a  dominant  role  in  the  overall 
water-induced  explosion-mitigation  process.  The  detonation -product  -to-water  momentum  transfer 
which  causes  water  aerosolization,  on  the  other  hand,  is  found  to  be  a  key  prerequisite  for  efficient 
explosion  mitigation. 
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1.  INTRODUCTION 

It  is  well  established  that  bulk  water  barriers  placed  in  the  vicinity  of  high-energy 
explosives  in  either  a  confined  environment  or  an  unconfined  environment  can 
significantly  mitigate  the  effects  of  explosions  by  reducing  the  quasi-static  gas  pressure 
and  impulse  (in  some  cases  up  to  90%)  [e.g.1-10].  This  explosion  mitigation  effect  of 
water  is  of  major  importance  for  the  safety  of  explosive  production  and  storage  facilities, 
protection  of  military  structures,  vehicles,  personnel  and  safe  demolition  of  employed 
unexploded  ordnance.  While  the  explosion  mitigation  effect  of  water  is  clearly 
demonstrated,  the  underlying  physical  mechanism  or  mechanisms  are  not  that  well 
understood.  For  example,  Keenan  and  Wager  [11]  suggest  that  blast  waves  from  the 
explosion  cause  the  water  to  aerosolize  and  upon  mixing  with  the  hot  detonation  gases  to 
evaporate.  This  causes  the  detonation  gases  to  cool  down  and  the  gas  pressure  to 
decrease.  In  addition,  the  resulting  lower  temperatures  of  the  gases  hamper  the  progress 
of  the  secondary  oxidation  reaction  of  the  detonation  products  in  oxygen- deficient 
explosives  further  reducing  the  gas  pressure.  Ericksson  and  Vretblad  [12],  on  the  other 
hand,  observed  that  the  latter  effect  may  not  be  significant  since  no  major  differences  in 
the  explosion  mitigation  effect  by  bulk  water  was  observed  in  highly  oxygen- deficient 
explosives  such  as  TNT  and  oxygen  non-deficient  explosives.  Also  a  frequently  cited 
mechanism  for  water-induced  explosion  mitigation  which  considers  bulk  water  as  a 
passive  barrier  that  can  absorb  some  of  the  energy  and  momentum  of  the  explosion 
products  is  not  well  understood  [13]. 

The  objective  of  the  present  work  is  to  carry  out  various  computational  analyses  in 
order  to  help  better  understand  and  quantify  the  contributions  of  various  mechanisms  to 
the  explosion  mitigation  effect  of  bulk  water.  It  should  be  noted  that  the  term  ‘ bulk 
water ”  is  used  throughout  the  paper  to  indicate  that  initially  water  is  present  in  a  bulk 
form,  however,  as  will  be  discussed  later,  the  interaction  of  the  explosion  products  with 
bulk  water  can  cause  the  water  to  aerosolize,  i.e.  break  up  into  micron-size  droplets. 
Also  the  term  ‘ bulk  water ”  is  used  to  make  a  distinction  relative  to  the  direct  use  of 
water  mist. 

The  organization  of  the  paper  is  as  follows.  A  brief  description  of  various 
mathematical  models  used  to  analyze  the  different  mechanical/thermal  and  chemical 
aspects  of  a  detonation  process  in  the  presence  of  bulk  water  and  an  overview  of  the 
constitutive  laws  governing  the  response  of  various  materials  when  subjected  to  blast 
loading  are  discussed  in  Sections  2. 1-2.3.  The  results  obtained  in  the  present  work  are 
presented  and  discussed  in  Section  3.  The  main  conclusions  resulting  from  the  present 
work  are  summarized  in  Section  4. 

2.  COMPUTATIONAL  ANALYSES 

2.1.  General  Consideration 

As  mentioned  earlier,  various  water-induced  explosion- mitigation  effects  are  analyzed 
computationally  in  the  present  work  using  several  analytical  and  numerical  techniques 
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ranging  from  those  emphasizing  thermo-chemical  aspects  of  the  detonation  process  to 
those  dealing  with  transient;  nonlinear-dynamics  based  mechanical  interactions  between 
detonation  products,  air  and  water.  In  general,  thermo -chemical  analyses  are  carried  out 
using  MATLAB,  a  general  purpose  mathematical,  computational  and  visualization 
package  [14].  Since  this  computer  package  is  widely  used  and  well-known,  it  is  not 
discussed  any  further  in  the  present  paper.  The  transient  nonlinear -dynamics  based 
mechanical  interactions  between  various  materials  participating  in  a  detonation  process 
are  analyzed  using  AUTODYN,  a  state-of-the-art  non-linear  dynamics  computational 
code  [15].  A  brief  overview  of  this  code  is  given  in  next  section. 

2.2.  Transient  Nonlinear- Dynamics  Calculations 

All  the  transient  nonlinear -dynamics  mechanical  interactions  between  various  materials 
attending  a  detonation  process  are  analyzed  computationally  using  AUTODYN,  a  state- 
of-the-art  nonlinear  dynamics  modeling  and  simulation  software  [15].  In  this  section,  a 
brief  overview  is  given  of  the  basic  features  of  AUTODYN,  emphasizing  the  aspects  of 
this  computer  program  which  pertain  to  the  problem  at  hand. 

AUTODYN  is  a  fully  integrated  engineering  analysis  computer  code  which  is 
particularly  suited  for  modeling  the  explosion,  blast,  impact  and  penetration  events. 
Codes  such  as  AUTODYN  are  commonly  referred  to  as  “ hydrocodes” .  Within  the  code, 
the  appropriate  mass,  momentum  and  energy  conservation  equations  coupled  with  the 
materials  modeling  equations  and  subjected  to  the  appropriate  initial  and  boundary 
conditions  are  solved.  The  numerical  methods  used  for  the  solution  of  these  equations 
involve  finite  difference,  finite  volume  and  finite  element  methods  and  the  choice  of  the 
method  used  (i.e.  ‘processor”  as  referred  to  in  AUTODYN)  depends  on  the  physical 
nature  of  the  problem  being  studied.  The  power  of  AUTODYN  is  derived  mainly  from 
its  ability  to  handle  complex  problems  in  which  different  regions  can  be  analyzed  using 
different  methods  such  as  the  Lagrange  processor  (typically  used  for  solid  continuum 
and  structures)  and  the  Euler  processor  (commonly  used  for  modeling  gases,  liquids  or 
solids  subject  to  large  deformations).  While  the  available  Euler  processor  provides 
multi-material  capabilities,  an  additional  Euler-FCT  single  material  processor  in  which 
materials  are  combined  to  a  single  material  using  a  Flux  Corrected  Transport  (FCT) 
approach  is  available  to  help  handle  computationally  intensive  multi- material  blast 
phenomena. 

Additional  methods  available  in  AUTODYN  include:  an  ALE  (Arbitrary  Lagrange 
Euler)  processor  capable  of  carrying  out  an  automatic  rezoning  (remeshing)  of  distorted 
grids;  a  Shell  processor  designated  for  modeling  thin  structures  and  a  gridless  SPH 
(Smooth  Particle  Hydrodynamics)  processor  which  does  not  suffer  from  a  grid  tangling 
problem  (typically  encountered  in  Lagrangian  processor)  and  does  not  entail  the  use  of 
an  unphysical  erosion  algorithm  (removal  of  highly  distorted  grids  to  help  the  numerical 
procedure). 

2.3.  Materials  Constitutive  Models 

Hydrodynamic  computer  programs  such  as  AUTODYN  [15]  are  capable  of  predicting 
an  unsteady,  dynamic  motion  of  a  material  system  by  solving  the  appropriate  mass, 
momentum  and  energy  conservation  equations,  subjected  to  the  associated  initial  and 
boundary  conditions.  However,  for  the  aforementioned  boundary  value  problem  to  be 
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fully  defined,  additional  relations  between  the  flow  variables  (pressure,  density,  energy, 
temperature,  etc.)  have  to  be  defined.  These  additional  relations  typically  involve  an 
equation  of  state,  a  strength  equation  and  a  failure  equation  for  each  constituent  material. 
These  equations  arise  from  the  hct  that,  in  general,  the  total  stress  tensor  can  be 
decomposed  into  a  sum  of  a  hydrostatic  stress  (pressure)  tensor  (which  causes  a  change 
in  the  volume/density  of  the  material)  and  a  deviatoric  stress  tensor  (which  is 
responsible  for  the  shape  change  of  the  material).  An  equation  of  state  then  is  used  to 
define  the  corresponding  functional  relationship  between  pressure,  density  and  internal 
energy  (temperature),  while  a  strength  relation  is  used  to  define  the  appropriate 
equivalent  plastic -strain,  equivalent  plastic  -strain  rate,  and  temperature  dependences  of 
the  equivalent  deviatoric  stress  (or  some  function  of  it).  In  addition,  a  material  model 
generally  includes  a  failure  criterion,  i.e.  an  equation  describing  the  (hydrostatic  or 
deviatoric)  stress  and/or  strain  condition  which,  when  attained,  causes  the  material  to 
fracture  and  loose  its  ability  to  support  normal  and  shear  stresses. 

In  the  present  work  the  following  materials  are  utilized  within  the  computational 
domain:  air,  a  high-energy  explosive  (e.g.,  TNT),  water  and  4340  steel.  In  the  following 
sections,  a  brief  description  is  given  of  the  models  used  for  each  of  the  four  constituent 
materials. 

2.3.1  Air 

Air  is  modeled  as  an  ideal  gas  and,  consequently,  its  equation  of  state  is  defined  by  the 
ideal-gas  gamma-law  relation  as  [15]: 

P  =  (y-l)-P-£  (1) 

Po 

where  P  is  the  pressure, y  the  constant-pressure  to  constant-volume  specific  heats  ratio 

(=1.4  for  a  diatomic  gas  like  air),  p0(=1.225kg/m  )  is  the  initial  air  density,  and  p  Is 
the  current  density.  For  Eq.(l)  to  yield  the  standard  atmosphere  pressure  of  101.3kPa, 
the  initial  specific  internal  energy  E  is  set  to  253.4kJ/m3  which  corresponds  to  the  air 
mass  specific  heat  of  717.6J/kg  K  and  a  reference  temperature  of  288.2K. 

Since  air  is  a  gaseous  material  and  has  no  ability  to  support  either  shear  stresses  or 
negative  pressures,  no  strength  or  failure  relations  are  required  for  this  material. 

2.3.2  High-Energy  Explosives 

High-energy  explosives  like  TNT  are  chemical  substances  which  when  subjected  to  a 
particular  mechanical  or  thermal  stimulus,  can  undergo  a  chemical  reaction  (or  a  series 
of  chemical  reactions)  releasing,  very  rapidly  (at  a  time  scale  of  microseconds),  a  large 
amount  of  energy.  The  whole  process  comprised  of  chemical  reactions  taking  place 
within  the  high-energy  explosive  and  the  associated  release  of  energy  is  generally 
referred  to  as  detonation.  Within  a  hydrodynamic  code,  detonation  is  typically  assumed 
to  take  place  instantaneously  at  each  material  point  within  the  high-energy  explosive 
when  such  points  are  properly  stimulated.  In  other  words,  as  a  result  of  initial  detonation 
at  predefined  detonation  locations  (points,  lines  or  planes),  detonation  waves  are 
generated  which  extend  into  the  unreacted  explosive  material  instantaneously  liberating 
the  energy  of  chemical  reactions  and  transforming  the  explosive  material  into  gaseous 
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detonation  products  at  each  point  swept  by  the  detonation  wave. 

Mathematical  description  of  a  detonation -induced  shock  wave  is  based  on  the  normal 
Rankine-Hugoniot  mass,  momentum  and  energy  conservation  equations  which  relate  the 
hydrodynamic  quantities  (density,  pressure,  energy)  across  the  discontinuity  (a  shock 
wave)  [16].  The  only  difference  between  the  Rankine-Hugoniot  equations  for  a  shock 
wave  in  a  non-explosive  material  and  their  counterparts  for  a  detonation  wave  is  the 
inclusion  of  a  chemical  energy  term  (initially  residing  as  a  chemical  energy  term  in  the 
explosive)  in  the  energy  equation. 

For  a  high-energy  explosive  material  in  an  initial  condition  defined  by  pressure  P0, 
density  p0  and  specific  volume  v0  =  l/p0,  the  line  denoted  as  ‘ Hugoniot  Curve”  in 

Fig.  1  represents  the  locus  of  all  (P,  v)  states  attainable  from  the  initial  (PQ,v0)  state 
during  detonation.  It  should  be  noted  that,  in  contrast  to  the  case  of  a  non -explosive 
material,  the  Hugoniot  curve  does  not  pass  through  the  initial  explosive  state  ( ,  V0 ) 
but  it  is  rather  shifted  toward  higher  pressures  due  to  the  chemical  energy  release 
accompanying  detonation  of  the  explosive  material.  In  addition  to  the  Hugoniot  curve, 
an  additional  line  denoted  as  the  ‘ Rayleigh  Line”  is  depicted  in  Fig.l.  This  line 
corresponds  to  the  all  (P,  \)  states  of  the  detonation  products  consistent  with  a  given 
detonation- wave  speed  and  the  Rankine-Hugoniot  mass  and  momentum  conservation 
equations.  It  should  be  noted  that  there  is  a  separate  Rayleigh  line  for  each  value  of  the 
detonation- wave  speed.  The  Rayleigh  line  depicted  in  Fig.  1  is  the  line  associated  with 
the  minimum  detonation  speed  which,  at  the  same  time,  yields  the  (P,  v)  states  of  the 
detonation- product  which  are  consistent  with  the  Hugoniot  curve.  The  point  of  tangency 
of  the  Rayleigh  line  and  the  Hugoniot  curve  is  generally  referred  to  as  the  Chapman- 
Jouget  point  and  defines  both  the  value  of  the  detonation  speed  and  the  corresponding  (P, 
v)  state  of  the  detonation  products. 


Fig.l.  Pressure  vs.  specific  volume  relations  in  a  high-energy  explosive  material. 
Please  see  the  text  for  details. 
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The  equation  of  state  for  an  explosive  then  defines  a  locus  of  (P,  v)  states  associated 
with  expansion  (or  compression)  of  the  detonation  products.  One  of  these  states  is 
clearly  the  initial  state  of  the  detonation  products  corresponding  to  the  Chapman-Jouget 
point.  Since  the  detonation  products  are  generally  associated  with  very  high 
temperatures,  their  expansion/compression  does  not  give  rise  to  a  significant  change  in 
entropy.  Hence,  it  is  generally  accepted  that,  at  least  in  the  vicinity  of  the  Chapman- 
Jouget  point,  the  equation  of  state  prescribes  the  states  which  are  consistent  with  the  P-v 
adiabat  passing  through  the  Chapman-Jouget  point.  The  remainder  of  the  P-v  curve 
corresponding  to  the  expansion  of  the  detonation  products  is  generally  obtained  by 
matching  the  area  under  the  P-v  curve  up  to  large  expansions  and  down  to  atmospheric 
pressure  with  the  experimentally  measured  work  of  detonation. 

The  Jones-Wilkins-Lee  (JWL)  equation  of  state  is  used  for  a  high-energy  explosive  in 
the  present  work  since  that  is  the  preferred  choice  for  the  equation  of  state  for  high- 
energy  explosives  in  most  hydrodynamic  calculations  involving  detonation.  The  JWL 
equation  of  state  is  defined  as  [17,  18]: 


p  =  A 


w 


-R.v 


>'  +Bt 
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-R,  i 


+ 


wE 
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where  the  constants  Ah  Rh  B h  R2  and  w  for  a  given  high-energy  explosive  like  TNT  are 
defined  in  the  AUTODYN  materials  library  and  v  is  the  specific  volume  of  the  material. 
As  explained  earlier,  within  a  typical  hydrodynamic  analysis,  detonation  is  modeled  as 
an  instantaneous  process  which  converts  unreacted  explosive  into  gaseous  detonation 
products  and  detonation  of  the  entire  high-energy  explosive  material  is  typically 
completed  at  the  very  beginning  of  a  given  simulation.  Consequently,  no  strength  and 
failure  models  are  typically  used  for  the  high-energy  explosives. 


2.3.3  Water 

Separate  equations  of  state  are  used  for  water  depending  on  whether  water  is  subjected 
to  expansion  or  compression.  When  water  is  subjected  to  expansion,  a  two-phase 
equation  of  state  proposed  by  Morgan  [19]  is  used  while  in  compression  a  polynominal 
type  of  equation  of  state  [15]  is  used. 

The  two-phase  equation  of  state  for  water  is  applicable  when  water  is  expanded  from 
its  initial  state  or  in  a  single -phase  (liquid  or  vapor)  region  or  in  the  two- phase 
(water+vapor)  region.  As  shown  in  Fig.2,  the  pressure/specific -volume  plane  at  the 
expansion  side  (  v  >  V0  ,  V0the  initial  specific  volume)  is  divided  into  two  regions  by  the 

saturation  curve.  Above  this  curve  water  exists  in  a  single -phase  (liquid  or  vapor)  state, 
while  below  the  saturation  curve  the  two  phases  coexist.  Values  for  P,  v,p  ,  e  and  T 

along  the  saturation  curve  are  known  and  are  tabulated  with  the  AUTODYN  material 
database. 

When  the  values  of  v  and  E  place  water  in  a  single -phase  region,  pressure  is  defined 
by  the  Mie-Gruneisen  equation  of  state  in  the  form  [15]: 

T  (v) 

P  =  P(v)  +  ^±l[E-Er(v)] 


v 


(3) 
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where  the  reference  quantities  denoted  with  a  subscript  r  are  known  and  are  associated 
with  the  saturation  curve,  while  F  =  v(dP/dE)v  is  the  so-called  Gruneisen  gamma, 
whose  values  are  also  known  along  the  saturation  line. 


Fig-2  .Pressure  vs.  specific  volume  relations  for  water  used  in  the  derivation  of  the  two- 
phase  equation  of  state.  Please  see  the  text  for  details. 


When  the  values  of  v  and  E  place  water  within  a  two-phase  region,  the  volume  and 
energy  weighted  average  equations: 

v  —  OCv  +  (1 — OC  )v i ,  and  (4) 

<5  1 


E  =  aE^  +(l-a  )E. 

<5  ' 


(5) 


where  a  is  the  mass  fraction  of  vapor  and  subscript  g  and  /  denote  the  gas  phase  and  the 
liquid  phase,  respectively,  are  combined  to  yield: 


v  —  v„ 


v„  -V, 


E~Es 

E„  -E, 


(6) 


Since  the  specific  volumes  and  energies  of  the  saturated  liquid  and  vapor  (v, ,  Eh 

Vg,  Eg  )  are  mutually  constrained  by  the  saturation  curve,  Eq.  (6),  defines  the  pressure  (or 

temperature)  for  the  given  values  of  v  and  E. 

When  water  is  subjected  to  compression,  the  following  polynomial  equation  of  state  is 
used: 


P  =a1\i  +  a2[i2  +a3|i3  +  (bQ  +bx\i)poE 


(7) 
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where  H  =  P/P0-1  is  the  compression,  p0  the  initial  density  and  the  coefficients  «/, 

a 2,  cit,  bo  and  bj  are  defined  in  the  AUTODYN  material  library. 

No  strength  model  is  used  for  water  because  water  has  very  little  ability  to  support  a 
shear  stress.  However,  in  order  to  examine  the  possibility  for  water  break-up  into  small 
droplets  as  a  result  of  its  interaction  with  the  detonation  gases,  a  minimum  negative- 
pressure  failure  criterion  is  used.  The  magnitude  of  the  critical  minimum  failure-pressure 
is  varied  in  order  to  examine  its  effect  on  the  results  of  the  computational  analyses. 


2.3.4  4340  Steel 


For  inert  solid  materials  like  4340  steel,  a  linear  type  of  equation  of  state  is  typically 
used  which  assumed  a  Hooke's  law  type  relationship  between  the  pressure,  P,  and  the 

r-  ^ 

volume  change  |4  ~  — — ~ 1  I  as: 


V 


J 


P=B\l 


(8) 


where  K  is  the  bulk  modulus  of  the  material.  Within  the  AUTODYN  material  database, 
the  initial  material  density  p0,  the  bulk  modulus  B.  the  specific  heat,  Cp  and  the 

reference  temperature  (Tref )  are  defined  for  4340  steel. 

To  represent  the  constitutive  response  of  4340  steel  under  deviatoric  stress,  he 
Johnson-Cook  model  is  used.  This  model  is  capable  of  representing  the  material 
behavior  displayed  under  large- strain,  high  deformation  rate,  high-temperature 
conditions,  of  the  type  encountered  in  problems  dealing  with  the  interactions  of 
detonation  products  and  solid  structures.  Within  the  Johnson-Cook  model,  the  yield 
stress  is  defined  as: 

r  =  k  +£,e;  Ji+c,bge  m 

where  8  is  the  equivalent  plastic  strain,  £p  the  equivalent  plastic  strain  rate,  Asr  the 
zero  plastic -strain,  unit  plastic -strain  rate,  room- temperature  yield  stress,  Bst  the  strain 
hardening  constant,  llst  the  strain  hardening  exponent,  Cst  the  strain  rate  constant,  mst 
the  thermal  softening  exponent  and  TH  =  (T  - Troom)l ( Tmell  - Tmm )  a  room  temperature 

( ^ room  )  based  homologous  temperature  while  Tmelt  is  the  melting  temperature.  All 
temperatures  are  given  in  degrees  of  Kelvin. 

Since  4340  steel  structures  are  generally  subjected  to  compressive  type  of  stresses,  no 
failure  model  was  used  for  the  4340  steel  in  the  present  work. 


3.  RESULTS  AND  DISCUSSION 


In  this  section,  the  four  most  frequently  cited  water -induced  explosion -mitigation 
mechanisms:  (a)  absorption  of  the  detonation  energy;  (b)  reduction  of  the  shock  speed; 
(c)  detonation-product-to-water  momentum  transfer  and  (d)  deceleration/suppression  of 
combustion  of  the  detonation  products  are  analyzed  using  the  underlying  physical 
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principles.  In  addition,  the  results  of  numerical  analyses  of  the  four  mechanisms  are 
presented  and  discussed.  In  each  case,  an  effort  was  made  to  quantify  the  effectiveness 
of  the  specific  water -mitigation  mechanism  in  reducing  the  pressure  (and  temperature) 
of  the  detonation-product  gas-mixture.  Such  effectiveness  is  judged  by  both  the 
magnitude  of  the  effect  and  by  the  rate  at  which  such  effect  takes  place. 


3.1  Absorption  of  the  Detonation  Energy 

3.1.1  General  Consideration 


Due  to  its  relatively  large  specific  heat  (~4.187kJ/kg  K)  and  a  relatively  large  latent  heat 
of  evaporation  (~2.25MJ/kg),  water  has  an  excellent  ability  to  absorb  energy  of 
detonation  (and  combustion)  of  high-energy  explosives.  For  comparison,  the  detonation 
energy  of  TNT  is  about  4.45MJ/kg  while  the  combined  energies  of  detonation  and 
subsequent  combustion  of  this  oxygen- deficient  explosive  s  about  4.45+10.22=14.67 
MJ/kg.  However,  for  the  water  to  be  able  to  absorb  the  detonation/combustion  energy  of 
a  high-energy  explosive  at  the  explosion  time  scale  of  tens  of  microseconds  to  tens  of 
milliseconds,  the  water  must  be  present  in  the  form  of  a  mist  (micron -size  droplets)  to 
ensure  a  fast  explosion- products -to-water  thermal  energy  transfer.  In  general,  the  mass 
evaporation  rate  of  water,  m  H  0  ,  can  be  defined  as: 


m 


H,0 


=  K 


A  (  psat  _  Y 
■H20  '  ‘h2o  a 


H.O 


P) 


(10) 


where  Km  is  the  liquid-to-vapor  mass  transfer  coefficient,  AH0  the  liquid/gas  surface 

area,  Pf0  the  water  vapor  saturation  pressure,  X  /A  0  the  molar  fraction  of  water- vapor 

in  the  gas  phase  and  P  the  pressure  of  the  gas  phase. 

According  to  Eq.  (10),  the  rate  of  evaporation  scales  linearly  with  the  surface  area  of 
water  being  evaporated.  The  surface  area  of  water  dispersed  into  N  spherical  droplets 
of  radius,  rH  0  ,  is  related  to  the  total  water  volume,  V H  0 ,  and  to  the  droplet  radius  as: 


A/,o  ~  ■ 


3  V, 
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The  application  of  Eq.  (11)  shows  that  a  single  water  sphere  which  has  a  volume  of  lm3 
and  a  surface  area  of  4.836nf  gives  rise  to  a  surface  area  of  -200, 000m  ,  if  it  is 
dispersed  into  spherical  droplets  with  a  radius  rH  O  =15(J,m.  In  accordance  with  Eq.  (10), 
this  would  increase  the  evaporation  rate  by  a  factor  of  over  41,000. 


3.1.2  A  Thermo-chemical  Analysis 

Evaporation  of  the  water  droplets  gives  rise  to  a  reduction  in  the  temperature,  and  hence, 
the  pressure  of  the  hot  detonation  products  and,  as  will  be  discussed  later,  can  slow 
down  or  even  prevent  the  subsequent  combustion  chemical  reaction  of  the  detonation 
products.  To  quantify  the  effect  of  water  evaporation  on  temperature  reduction  of  the 
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detonation  products,  a  simple  quantitative  thermal  energy  balance  analysis  for  the  case 
of  TNT,  a  high-energy  explosive  is  carried  out  in  this  section. 

TNT  is  a  solid  high-energy  explosive  with  a  chemical  formula  C7H5N3O6  which 
produces  during  detonation  a  gas  mixture  which  is  laden  with  solid-carbon  soot  and 
which  has  the  following  chemical  composition  in  mole  fractions:  X H  0  =  0.227  , 

X co  =0.318  ,  X C(sj  =0.318  and  X  Vi  =0.136  .  The  initial  temperature  of  the 

mixture  is  about  3500K.  Temperature  variation  of  the  molar  specific  heats  for  the  gas 
mixture  constituents  (as  well  as  for  molecular  oxygen  and  carbon  dioxide)  are  displayed 

in  Fig. 3.  From  the  molecular  weights  of  the  mixture  constituents,  M H  0  =  18  , 
Mco  =  28  ,  A7C(  v|  =12  and  MN  =  28  ,  all  in  g/mol,  the  average  molar  specific  heat, 
Cp,  u,  mix,  and  the  average  molecular  weight,  Mmj  ,  of  the  mixture  are  computed  as  a 
mole-fraction  weighted  averages  of  the  respective  constituent  quantities  as: 

C^,„=Ex,C,„,  (i  =  H20,C0, C(s), N2 )  (12) 

i 

and 

Mnax  =  Y<XiMi  ( i  =  H20,C0,C(s),N2 )  (13) 


The  corresponding  mass-based  specific  heat  of  the  mixture  is  then  obtained  as: 
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P  ,m,  mix 
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(14) 


To  quantify  the  effect  of  water  evaporation  on  temperature  reduction  of  the 
detonation- product  gas  mixture,  the  following  energy  balance  equation  is  used: 

T final 

mH 2o  '  CPm  H2o  (373  —  298)  +  mH  0  ■  Qevap,H2o  =  J '  m,,,ix '  CPjn  mixdT  (15) 

3500 


The  left-hand  side  in  Eq.  (15)  corresponds  to  the  energy  absorbed  by  a  m//  0  mass  of 

water  from  its  initial  temperature  (298K)  till  the  point  of  boiling  (T=373K)  and  during 
subsequent  evaporation.  The  right-hand  side  of  Eq.  (15),  on  the  other  hand,  corresponds 

to  the  energy  given  off  by  a  Mmjx  mass  of  the  hot  detonation -product  gas-mixture  and  is 
defined  to  take  into  account  temperature  dependence  of  the  mixture  specific  heat, 
mix,  in  accordance  with  Fig. 3. 

When  Eqs.  ( 12)  -  ( 15)  are  applied  to  one  kilogram  of  TNT  and  one  kilogram  of  water 
and  the  following  data  are  used  for  the  mass  specific  heat  and  the  mass  evaporation 

energy  of  water:  CP  m  H^0  =  4. 187 kJ /  kg  ■  K  and  QevapMl0  =  2.25 MJ  /  kg  ,  the  final 

temperature  of  the  mixture  is  obtained  asTfinal  =  2298 K  .  This  finding  shows  that,  if 
one  kilogram  of  water  can  be  evaporated  for  each  kilogram  of  TNT  detonated,  the 
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temperature  of  the  resulting  detonation  gas  mixture  would  be  reduced  by  ?  7=3500  - 
2298  =1202K. 


Fig.3.  Temperature  variation  of  the  molar  specific  heat  for  a  number 
of  species  attending  detonation/combustion  of  TNT. 

Clearly,  water  can  be  very  efficient  in  reducing  the  temperature  of  the  detonation- 
product  gas -mixture.  In  fact,  since  the  specific  heat  of  the  gas -mixture  constituents 
decreases  as  their  temperature  is  reduced  in  Fig.3,  further  increase  in  the  relative  amount 
of  water  would  be  even  more  effective  in  reducing  the  temperature  of  the  detonation- 
product  gas -mixture.  As  will  be  discussed  later,  a  reduction  in  the  temperature  of  the 
detonation- product  gas -mixture  can  substantially  reduce  the  rate  of  subsequent  oxidation 
(combustion)  of  CO(g)  and  C(s)  in  the  mixture  or  even  prevent  this  reaction  from  taking 
place.  Since,  this  reaction  is  associated  with  the  formation  of  a  fire-ball  which  advances 
with  the  detonation  front  and  releases  a  substantial  amount  of  thermal  energy,  its 
suppression  is  highly  desirable  from  the  explosion-mitigation  stand  point. 

Previous  experimental  investigations  [e.g.  HO]  have  shown  that  the  use  of  water 
barriers  can  reduce  the  peak  hydrostatic  pressure  by  as  much  as  90%.  The  computed 
water-induced  reduction  in  the  temperature  of  the  detonation- product  gas-mixture  from 
-3500K  by  -1200K,  would  certainly  give  rise  to  a  reduction  n  the  pressure  of  the 
mixture.  However,  it  is  not  clear  that  the  reduction  in  the  mixture  temperature  could 
solely  account  for  the  observed  reduction  of  the  peak  hydrostatic  pressure.  If  one 
assumes  that  the  hydrostatic  pressure  of  a  gas  mixture  scales  linearly  with  its 
temperature,  as  is  the  case  for  the  ideal  gas  law  (the  use  of  ideal  gas  law  might  not  be 
fully  appropriate  considering  high  pressure  levels  in  the  detonation -product  gas-mixture) , 
the  observed  temperature  reduction  would  give  rise  to  a  pressure  reduction  in  a  range 
between  30-35%.  While  increasing  the  relative  amount  of  water  is  expected  to  give  rise 
to  a  further  reduction  in  the  gas-mixture  temperature  and  pressure,  the  extent  of  this 
reduction  may  not  be  proportional  to  the  amount  of  water  since  it  depends  on  the  ability 
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of  the  detonation-product  gas-mixture  to  aerosolize  the  additional  water. 

As  stated  above,  the  use  of  the  ideal  gas  law  to  predict  the  potential  water- evaporation 
induced  reduction  in  the  gas -phase  pressure  may  not  be  very  reliable  due  to  the  attendent 
high  pressure  levels.  Here  we  provide  yet  another  estimate  of  this  pressure  reduction. 
In  several  experimental  investigations  [e.g.  10],  it  was  established  that  the  peak  side-on 
(hydrostatic)  pressure  at  a  given  point  scales  with  the  explosive  total  energy  raised  to  a 
power  of  1/3.  If  one  treats  the  energy  absorbed  by  water  during  evaporation  as  a  means 
of  reducing  the  total  energy  of  the  explosive,  one  can  obtain  another  estimate  for  the 
water-induced  pressure  reduction. 

As  stated  earlier  the  total  energy  of  explosion  for  TNT  which  includes  the  heat  of 
detonation  and  the  heat  of  subsequent  combustion  is  ~  14.67MJ/kg  TNT.  Evaporation  of 
lkg  of  water  consumes  about  2.25MJ  of  energy,  so  that  the  effective  total  energy  of  TNT 
can  be  considered  as  14.67-2.25=1 2.42MJ/kg  TNT,  a  reduction  of  -15%  which  would 
yield,  according  to  the  aforementioned  correlation,  a  decrease  in  pressure  of  only  ~  5% 
per  kg  of  water.  If  water  completely  suppresses  the  combustion  reaction,  on  the  other 
hand,  the  effective  specific  energy  of  TNT  is  reduced  to  14.67-10.22-2.25=2.2MJ/kg 
TNT,  a  reduction  of  -  85%.  The  resulting  pressure  reduction  is  about  53%.  It  appears, 
hence,  that  for  the  water  evaporation  process  to  have  a  significant  effect  on  gas-phase 
pressure  reduction,  it  must  cause  cessation  of  the  combustion  reaction  in  oxygen- 
deficient  explosives. 

It  should  be  noted  that  aerosolization  of  water  also  consumes  some  energy  initially 
residing  in  the  detonation -product  gas -mixture.  Namely,  the  formation  of  small  water 
droplets  is  accompanied  by  a  substantial  increase  in  the  energy  stored  at  the  surface  of 
water.  The  surface  energy  of  water  is  0.072J/nT  at  room  temperature  and  decreases  to 
~0.060J/nr  near  the  boiling  point.  As  shown  in  the  previous  section,  lm3  of  water 
dispersed  into  1 5  pm- radius  droplets  have  a  total  surface  area  of  ~  200,000m2.  Thus  the 
energy  absorbed  by  water  in  the  aerosolization  process  is  at  most  0.072x  200,000 
=14,400J/kg  FLO  =  14.4  kJ/kg  HO.  This  value  is  only  a  tiny  fraction  of  the  water 
evaporation  energy  (~2.25MJ/kg  H20)  and,  hence,  can  be  neglected. 

3.1.3  A  Fluid  -dynamic  Analysis 

According  to  Eq.  (10),  the  evaporation  rate  of  water  is  proportional  to  a  temperature- 
dependent  liquid -to- vapor  mass  transfer  coefficient,  Km.  To  a  great  extent,  the  magnitude 
of  Km  is  controlled  by  the  rate  at  which  the  heat  is  transferred  from  the  detonation- 
product  gas-mixture  to  the  water  droplets.  The  gas-mixture-to-water  droplets  heat 
transfer  rate  is,  on  the  other  hand,  greatly  affected  by  the  extent  of  turbulence  in  the  gas- 
phase  stream.  Furthermore,  the  presence  of  water  droplets  in  the  gas -phase  stream  is 
expected  to  have  some  effect  on  the  extent  of  turbulence.  In  the  early  stages  of  the  water 
break-up  process,  when  large  water  droplets  (fragments)  are  present,  a  general  increase 
in  the  extent  of  turbulence  is  expected  since  water  droplets  can  create  turbulent  regions 
in  their  wake.  Conversely,  in  the  later  stages  of  the  water  break-up  process,  when 
micron-size  water  droplets  are  present,  such  droplets,  due  to  their  large  density,  would 
reduce  flow  fluctuations  and,  hence,  lower  the  extent  of  turbulence.  However,  one  must 
take  into  account  the  fact  that  due  to  high  values  of  the  evaporation  rates,  the  life  span  of 
micron-size  droplets  is  very  short.  While  no  quantitative  analysis  is  carried  out  in  the 
present  work  pertaining  to  the  effects  of  water-droplets  number  density  and  size  on  the 
gas-to-water  heat  transfer  coefficient,  the  discussion  given  earlier  in  this  section 
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(specifically  the  presence  of  counter  balancing  phenomena)  suggests  that  their  effects 
are,  perhaps,  of  a  second  order  relative  to  the  effect  of  enhanced  surface  area  on  the  rate 
of  evaporation  of  the  water. 

3.1.4  A  Kinetic  Analysis  of  Water  Evaporation 

The  analyses  presented  above  demonstrate  that  water  evaporation  has  the  potential  for 
absorbing  the  substantial  portions  of  the  thermal  energy  carried  by  the  detonation- 
product  gases.  However,  for  this  potential  to  be  fully  utilized,  water  must  be  dispersed 
in  the  form  of  small- size  droplets  which  can  evaporate  over  a  time-period  of  up  to  few 
tens  of  microseconds.  In  this  section,  we  utilized  several  existing  droplet  evaporation 
models  to  assess  the  time  for  evaporation  of  the  water  droplets  of  different  sizes 
dispersed  within  a  detonation- product  gas  mixture  at  a  temperature  of  3500K. 

The  first  model  utilized  was  proposed  by  Butz  et  al.[26],  and  assumes  that  the  droplets 
do  not  move  relative  to  the  gas,  that  the  droplet  surface  temperature,  and  the  gas  far¬ 
away  temperatures  are  constant.  These  conditions  correspond  to  the  Nusselt  number  of 
2.  The  results  pertaining  to  the  droplet  diameter  dependence  of  the  evaporation  time,  are 
displayed  in  Fig. 4,  the  curve  labeled  "‘Butz  et  al.[26]”. 


Fig-4.  Variation  of  the  simple  water-droplet  evaporation  time 
with  the  droplet  diameter.  Please  see  text  for  detail. 

The  second  model  was  developed  in  the  present  work  and  is  briefly  described  in 
Appendix  A.  This  model  is  an  extension  of  the  model  developed  by  Butz  et  al.  [26],  and 
includes  the  effect  of  droplet  heating  to  the  boiling  temperature.  The  results  of  this 
analysis  are  also  displayed  in  Fig. 4,  the  curve  labeled,  “ Present  Work” . 

The  last  model  utilized  is  in  fact  a  series  of  non -equilibrium  evaporation  models 
described  by  Miller  et  al.  [25].  A  brief  description  of  this  model  is  given  in  Appendix  B. 
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The  predictions  given  by  this  model  are  shown  in  Fig. 4,  the  band  labeled  ‘ Miller  et  al. 

[25 r. 

The  results  displayed  in  Fig.  4,  suggest  that  in  order  for  water  droplets  to  be  able  to 
evaporate  within  a  time  span  not  greater  than  10ms,  their  diameter  should  be  smaller 
than  ~90|J,m.  The  period  of  10ms  is  commonly  considered  as  a  time  span  over  which  the 
major  portion  of  the  momentum  transfer  from  the  detonation -product  gas -mixture  to  the 
target  structure/personal  takes  place. 

3.2  Reduction  in  the  Shock  Speed 
3.2.1  General  Consideration 

When  water  is  aerosolized  as  a  result  of  an  interaction  with  the  detonation  products  of  a 
nearby  high-energy  explosive,  the  speed  of  sound  and,  thus,  the  shock  speed  in  the 
surrounding  moisture- laden  air  becomes  substantially  reduced  [20].  This  observation  can 
be  rationalized  as  follows: 

(a)  The  speed  of  sound  in  a  material  can  be  defined  as  the  square  root  of  a  ratio  of  its 
bulk  modulus,  B  ,  and  its  density,  p  ; 

(b)  At  room  temperature,  water  has  a  relatively  large  value  of  the  bulk  modulus  of 
Bh  o  =  2.05 X 1 0  ’  Pa  and  a  density  of  approximately  p H^Q  =  999kg / m3 ,  and 
in  accordance  with  the  aforementioned  relation,  the  speed  of  sound 
is  CH  0  =  1,432 ml  s\ 

(c)  In  the  frequency  range  of  the  sound  waves  and  the  shock  waves  accompanying  a 
detonation  event,  the  bulk  modulus  of  a  gas  is  typically  defined  as  the  product  of  its 
pressure,  P  ,  and  the  ratio  of  its  constant-pressure  and  its  constant-volume  specific 
heats,  y  .  For  a  diatomic  gas  like  air,  y=1.4,  and  at  the  atmospheric  pressure 

P=101.3kPa  and  at  a  standard  density  ofpni>  =1.225  kg/m3,  the  bulk  modulus 

and  the  speed  of  sound  are:  Bair  —  1 .414  X 10 5  Pa  and  Cair  —  340 .3m/ s  , 
respectively; 

(d)  When  water  is  finely  dispersed  in  air,  the  speed  of  sound  as  a  function  of  the 
volume  fraction  of  water  shows  a  very  interesting  behavior  as  displayed  in  Fig. 5. 
This  behavior  arises  from  the  fact  that  the  bulk  modulus  of  the  air -water  mixture, 

is  given  as: 


1 

B 

mix 

where  fHM  is  the  volume  fraction  of  water  in  the  mixture.  Eq.  (16)  reflects  the  fact  that 

the  modulus  of  the  air-water  mixture  is  dominated  by  the  phase  present  in  a  larger 
volume  fraction  and  by  the  phase  that  is  volumetrically  more  compliant. 

The  density  of  the  air -water  mixture,  pmix,  is  likewise  defined  as: 


1  f h2o  f h.o 

-  +  -  - 


B 

air 


B 


(16) 


H,0 


P  mix  =  (l-/«2o)P  air  f  h2o  P  h2o 


(17) 
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The  variation  of  the  speed  of  sound  in  an  air-water  mixture,  C  ,  with  the  volume 
fraction  of  water,  fHl<0  ,  for  the  previously  stated  values  of  B H  0 ,  Bair ,  p  H  0  and  pair, 
is  calculated  using  the  equation  Cmix(  fH2o  )  =  -jBmi  J  p  mix  . 

(e)  The  results  of  this  calculation  are  shown  in  Fig.5.  The  results  displayed  in  Fig.5 
show  that  the  speed  of  sound  in  an  air  -water  mixture  is  substantially  reduced  and  it 
is  below  50m/s  for  the  water  volume  fraction  between  0.06  and  0.94;  and 

(f)  According  to  the  standard  normal  shock  relations,  a  fixed  pressure  discontinuity 
across  a  shock  wave,  P2/P1,  corresponds  to  a  unique  value  of  the  Mach  number  for 
air  or  air-water  mixture  behind  the  shock  wave,  Ms.  A  fixed  value  of  the  Mach 
number  at  a  reduced  level  of  the  speed  of  sound  implies  a  lower  speed  of  air- water 
mixture  moving  behind  the  shock  wave  as  well  as  a  lower  speed  of  the  shock  wave. 


Fig.5.  The  variation  of  sound  speed  in  the  air/water  mixture  as  a  function 
of  the  volume  fraction  of  water  at  the  ambient  temperature. 


3.2.2  A  Nonlinear-dynamics  Shock-based  Analysis 

To  demonstrate  computationally  that  the  presence  of  moisture  in  air  would  indeed 
reduce  the  shock  speed,  a  simple  Euler -based  simulation  of  a  shock- tube  experiment  is 
carried  out  using  AUTODYN.  Towards  that  end,  a  two  dimensional  Eulerian 
axisymmetric  computational  domain  with  a  length-to-diameter  ratio  of  20:1  is  used. 
One  twentieth  of  the  tube  length  at  one  of  its  ends  is  tilled  with  high-pressure  air  while 
the  remainder  of  the  tube  is  tilled  with  either  the  atmospheric  -pressure  air  or  with  an  air- 
water  mixture  containing  50  volume  %  of  water  at  the  atmospheric  pressure.  Along  all 
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the  edges  of  the  computational  domain,  the  “ no-flow ”  boundary  conditions  are  applied. 
Six  gage  points  at  equal  intervals  of  50mm  were  placed  and  the  pressure  at  these  points 
recorded  during  the  simulation. 

The  results  pertaining  to  the  pressure-time  traces  at  the  six  gage  points  for  air  and  air- 
water  mixture  are  displayed  in  Figs  6(a)  -  (b),  respectively.  From  the  shockwave  arrival 
times  at  the  different  gage  points  (the  times  at  which  the  overpressure  reaches  the  peak 
value)  displayed  in  Figs  6(a) -(b)  and  the  known  locations  of  the  gage  points,  the  average 
shockwave  speeds  (and  their  temporal  evolutions)  are  computed.  The  ratio  of  the  two 
average  speeds  is  found  to  vary  between  14.0  aid  14.1  as  a  function  of  time.  A  similar 
ratio  is  found  for  the  longitudinal  velocities  of  air  and  the  air -water  mixture  behind  the 
shock  wave.  These  are  in  excellent  agreement  with  the  corresponding  ratio  of  the  sound 
speeds  obtained  using  the  procedure  described  in  the  previous  section 

Cair/Cniix  =  340.3m/ 5 /23.82m/  s  -14.3  . 

The  results  obtained  clearly  show  that  the  speed  at  which  the  air -water  mixture  moves 
and  ultimately  interacts  with  a  structure/personal  is  reduced.  This  speed  reduction,  on 
one  hand,  yields  longer  times  over  which  water  break  up  into  droplets  and  evaporation 
can  take  place.  On  the  other  hand,  as  demonstrated  in  Section  3.3.2,  this  also  causes  a 
reduction  in  the  relative  droplet/gas  velocity,  the  parameter  which  greatly  affects  the  rate 
of  the  water  break-up  process.  The  presence  of  these  two  counteracting  effects  suggests 
(in  the  absence  of  a  detailed  numerical  analysis)  that  the  shock-speed  reduction  has, 
perhaps,  a  second-order  effect  on  the  water  break-up  and  evaporation  processes . 

The  results  displayed  in  Figs  6(a)- (b),  show  that  the  presence  of  water  increases  the 
magnitude  of  the  hydrostatic  overpressure  by  around  20%.  Even  more  serious  is  the 
potential  increase  in  the  dynamic  pressure  caused  by  the  use  of  water.  The  dynamic 
pressure  generally  scales  with  a  product  of  the  density  of  the  fluid  and  the  square  of  its 
speed.  In  the  present  case,  the  use  of  50  volume  %  of  water  causes  the  density  to 
increase  by  a  factor  of  nearly  500  while  the  squared  velocity  would  decrease  by  a  factor 
of  14.3“=225.  This  finding  suggests  that  the  dynamic  pressure  would  in  fact  increase  by 
a  factor  (=500/225)  greater  than  2.  It  appears,  hence,  that  for  the  water  to  display  the 
desired  explosion- mitigation  effect,  it  should  be  dispersed  into  fine  droplets  but  these 
droplets  should  evaporate  (before  reaching  the  target)  and  reduce  the  temperature  and 
side -on  pressure  of  the  detonation- product  gas  -mixture  rather  than  remain  as  liquid 
droplets  dispersed  within  the  gas  phase. 

To  demonstrate  directly  that  water-aerosolization  induced  decrease  in  the  shockwave 
speed  may  not  have  the  often-predicted  explosion-mitigation  effect,  a  small  cylindrical 
disk  shaped  target  made  of  4340  steel  is  inserted  into  the  shock  tube.  The  axis  of  the 
target  is  set  to  coincide  with  that  of  the  tube.  The  target  is  modeled  as  a  Lagrangian  part 
and  the  Euler/Lagrange  coupling  option  inside  AUTODYN  is  used  to  account  for  the 
interaction  between  the  detonation-product  gas -mixture  (with  and  without  water  droplets 
dispersion)  and  the  target.  The  temporal  evolution  of  the  momentum  of  the  target  is 
monitored  in  the  two  cases.  It  is  found  that  the  target  momentum  is  approximately  twice 
as  large  in  the  case  of  the  explosion- products  gas-mixture  containing  water  droplets,  in 
full  agreement  with  the  results  of  the  dynamic -pressure  analysis  presented  above. 
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3.3  Momentum  Transfer 

3.3.1  General  Consideration 

When  the  detonation  front  of  a  high-energy  explosive,  explosive  casing  fragments  or 
soil  ejecta  reach  the  explosive/water -barrier  interface,  they  creates  a  positive  pressure 
wave  which  travels  through  the  water.  When  this  wave  reaches  the  opposite  water 
surface,  i.e.  the  water/air  interface,  it  reflects  as  a  negative-pressure  (expansion)  wave 
due  to  a  large  acoustic  impedance  mismatch  between  water  and  air.  Acoustic  impedance 
is  typically  defined  as  a  product  of  the  material  density  and  its  sound  speed  and  based  on 
the  density  and  sound  speed  values  for  air  and  water  reported  in  Section  3.2.1,  the 
acoustic  impedances  of  water  and  air  are  1.43xl()6Rayles  and  417.7Rayles,  respectively. 
Thus,  there  is  a  significant  acoustic  impedance  mismatch  at  the  water/air  interface. 

A  liquid  like  water  cannot  sustain  a  significant  negative  pressure  and  the  water  region 
at  the  water/air  interface  breaks  up  into  droplets  and  sprays  out  into  the  surrounding  air. 
This  process  continues  as  the  expansion  wave  propagates  back  towards  the 
explosive/water  interface.  In  the  process  just  described,  the  momentum  initially  carried 
by  the  explosive  products,  soil  ejecta  or  other  debris  are  transferred  to  the  effectively 
entire  mass  of  the  surrounding  water.  Consequently,  the  average  velocity  at  which  the 
water  droplets  are  sprayed  out  can  be  a  very  small  fraction  of  the  shock  velocity, 
substantially  mitigating  the  effects  of  explosion  on  the  surrounding  structures/personnel. 

3.3.2  Water  Fragmentation  Computational  Analysis 

To  demonstrate  computationally,  that  the  water  aerosolization  process  described  above  is 
feasible,  a  simple  SPH  (Smooth  Particle  Hydrodynamics)  method  based  model  is 
developed  and  solved  using  AUTODYN.  The  SPH  solver  is  used  because  its  non-mesh 
(i.e.  particle)  based  nature  gives  rise  to  a  more  clear  visual  description  of  the  water 
aerosolization  process. 

The  rectangular  computational  domain  used  in  this  portion  of  the  work  is  shown  in 
Fig. 7  (a).  Within  this  domain,  a  semicircular  region  is  tilled  with  TNT  while  the 
remainder  of  the  domain  is  filled  with  water.  The  lower  edge  of  the  domain  is  defined  as 
the  axis  of  symmetry,  while  “ zero-stresC  boundary  conditions  are  applied  to  the 
remaining  domain  edges.  A  single  detonation  point  is  placed  at  the  center  of  the  bottom 
edge.  Two  gage  points  are  placed  along  the  vertical  midline,  as  shown  in  Fig.7  (a)  so  that 
a  temporal  variation  of  the  pressure  can  be  monitored  at  these  points.  In  addition,  field 
variations  of  the  pressure  and  the  morphology  of  the  computational  domain  during 
simulation  are  also  monitored. 

Typical  results  pertaining  to  the  field  and  temporal  variation  of  the  pressure  show  a 
circular  (compressive)  shock  wave  emanating  from  the  TNT  detonation  point  and 
propagating  toward  the  opposite  long  edge  of  the  domain,  i.e.  the  water/air  interface. 
Once  the  shock  wave  reaches  the  opposite  long  edge  of  the  domain,  it  reflects  as  an 
expansion  wave.  This  behavior  is  more  clearly  seen  by  monitoring  the  pressure-time 
traces  a  the  two  gage  points,  as  shown  in  Fig.7  (b).  The  results  shown  in  Fig.7(b)  show 
that  the  shock  wave  arrives  first  at  the  gage  point  1  which  is  closer  to  the  TNT  charge, 
while  the  expansion  wave  first  reaches  the  gage  point  2  which  is  closer  to  the  water/air 
interface. 

The  results  pertaining  to  the  morphology  of  the  computational  domain  are  found  to  be 
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fairly  sensitive  to  the  negative -pressure  level  selected  for  the  onset  of  failure  in  water.  In 
general  however,  two  different  aerosolization  mechanisms  are  observed:  (a)  a  “ spalling' ’ 
mechanism  in  which  distinct  layers  of  water  is  spalled-off  at  the  water/air  interface.  Fig. 
8(a)  and;  (b)  a  “ water-fragmentation ”  process  in  which  water  at  the  water/air  interface 
breaks  up  into  irregular  fragments.  Fig. 8  (b).  The  former  mechanism  is  typically 
observed  when  the  least  negative  values  of  the  failure-pressure  are  used  while  the  latter 
is  typically  observed  at  the  most  negative  levels  of  the  failure  pressure.  In  the  case  of 
water  spallation,  discrete  layers  of  water  are  formed  and  the  interlayer  gaps,  initially 
containing  vacuum,  are  quickly  filled  with  water  vapor.  This,  in  turn,  prevents  water 
layers  from  coming  back  into  a  direct  contact  and  further  facilitates  break-up  of  the 
layers  into  discrete  water  droplets.  In  the  case  of  the  water-fragmentation  process,  water 

fragments  freely  fly  out  and  undergo  further  break-up  into  small-size  droplets. 

I 
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Fig.7.  (a)  A  computational  domain  used  in  the  SPH  (smooth  Hydrodynamics)  analysis  of  the 
water  aerosolization  process;  (b)  Typical  variation  of  pressure  at  the 
two  gauge  points  during  simulation. 
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Fig-8.  Two  mechanisms  of  water  aerosolization:  (a)  a  “spallation”  mechanism  and 
(b)  a  “fragmentation”  mechanism. 

The  two  mechanisms  described  above  can  be  considered  as  the  mechanism 
controlling  the  initial  stage  of  water  break-up  into  larger  fragments.  Once  these 
fragments  are  formed,  they  undergo  further  break-up  into  smaller  water  droplets.  A 
detailed  analysis  of  the  subsequent  fragment  break-up  process  is  not  carried  out  in  the 
present  work.  Nevertheless,  it  is  well  established  that  there  is  a  number  of  such  fragment 
break-up  mechanisms  including  bag  break-up,  multi-mode  break-up  and  shear  stripping 
[21].  A  simple  schematic  of  these  three  fragment  break-up  mechanisms  is  displayed  in 
Figs  9(a) -(c).  In  each  case,  a  spherical  droplet  (far  left)  is  deformed  as  a  result  of  its 
interaction  with  the  aerodynamic  forces  of  the  surrounding  gas.  In  the  case  of  the  bag 
break-up  mechanism,  deformed  droplets  acquire  a  thin-wall  half-ellipsoid  bag  shape. 
During  the  multi-mode  break-up  multiple  bags  are  formed  and  connected  to  a  central 
column  of  fluid  forming  an  umbrella  -  like  shape.  Within  the  shear -stripping  mode,  only 
the  outer  layers  of  the  droplets  are  deformed  (sheared  off).  Which  of  these  mechanisms 
dominates  water  breakup  into  small  droplets  is  controlled  by  the  magnitude  of  the  Weber 
number.  We,  which  is  defined  as: 


G 


(18) 


where  pw>  is  the  water  density,  ?  the  velocity  of  the  gas  stream  relative  to  that  of  the 
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water  fragment,  d  the  equivalent  droplet  diameter  and  s  the  water/gas  surface  energy.  It 
is  generally  recognized  [21-22]  that  the  minimum  critical  value  for  the  Weber  number 
required  for  the  break-up  mechanism  to  become  operative  is  about  12.  Multi-mode 
break-up  and  shear  stripping  mechanisms  entail  high  Weber  numbers  of  approximately 
20  and  70,  respectively. 

Since  the  water  density,  pw,  and  its  surface  tension  can  vary  only  in  relatively  small 
ranges,  large  values  of  the  Weber  number  are  attained  according  to  Eq.(18)  at  large 
values  of  the  relative  velocity  of  the  gas  and  at  large  values  of  the  equivalent  droplet 
diameter.  This  indicates  that  the  efficiency  of  the  water  break-up  process  is  the  highest  at 
the  shortest  post-detonation  times  and  degrades  continuously.  The  attainment  of  a  fine 
water  mist  as  a  result  of  an  interaction  between  the  detonation- product  gas-mixture  and 
the  bulk  water  is  governed  by  the  magnitude  of  the  explosive’s  internal  energy,  since  the 
explosives  internal  energy  controls  the  initial  detonation -product  gas  velocity.  Once  a 
water  mist  consisting  of  micron -size  droplets  is  produced,  as  argued  in  Section  3.2,  the 
evaporation  rate  can  become  significant,  causing  a  major  absorption  of  the  energy 
carried  by  the  detonation-product  gas-mixture  and  a  substantial  reduction  in  the  gas- 
phase  temperature  and  pressure.  In  other  words,  the  momentum  transfer  from  the 
detonation- product  gas -mixture  to  the  bulk  water  and  the  subsequent  aersolization  of  the 
water  are  the  key  prerequisites  for  obtaining  the  desired  explosion-mitigation  effect. 
The  actual  explosion -mitigation  effect,  however,  appears  to  be  dominated  by  the  water- 
droplets  evaporation  process.  It  is,  hence,  critical  that  the  water  break-up  process  into 
micron-size  droplets  is  completed  over  a  short  time  period  (perhaps  less  than  1ms)  to 
enable  a  sufficient  time  for  droplet  evaporation.  Based  on  the  substantial  amount  of 

experimental  data,  the  following  correlation  for  the  time,  tb  ,  for  one  break-up  event  has 
been  proposed  [27]: 


t 


b 
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(19) 


where  p  is  the  gas  density. 

As  established  earlier,  the  Weber  number  continuously  increases  as  a  function  of  time 
following  detonation,  so  that  the  time  for  each  subsequent  break-up  event  becomes 
longer  and  longer.  The  increase  in  the  Weber  number  during  simulation  is  being  studied 
by  us  in  an  ongoing  investigation.  Based  on  the  preliminary  results  of  this  investigation, 
and  under  the  assumption  that  the  initial  water  fragments  have  an  effective  diameter  of 
10mm,  the  time  required  to  produce  water  mist  with  an  average  droplet  diameter  of 
5 On m  has  been  estimated  between  0.5  and  2ms.  While  more  work  is  needed  in  the 
ongoing  analysis  to  narrow  down  this  time  range,  the  initial  results  appear  encouraging, 
with  respect  to  obtaining  micron-size  water  droplets  over  a  time  period  which  is 
relatively  small  compared  with  the  droplet  evaporation  time. 
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(c) 

Fig.9.A  schematic  of  three  water  break-up  mechanisms  commonly  observed: 
(a)  Bag  break-up;  (b)  Mixed-mode  break-up;  and  (c)  Shear  stripping. 


3.4  Deceleration/Suppression  of  Combustion 
3.4.1  General  Consideration 

Many  explosives  like  TNT  are  oxygen  deficient,  i.e.  they  do  not  contain  enough  oxygen 
to  oxidize  all  the  elements  in  the  explosive  during  an  explosion.  The  detonation  reaction 
of  TNT  with  a  chemical  formula  QH5N306  in  the  absence  of  external  oxygen  can  be 
written  as  follows; 

C7H5N306(s)  -+2.5H 20(g)  +3.5C(s)  +  3.5CO(g)  +  L5N2(g)  +1.034MJ / moleTNT {nQ) 

If  at  least  5.25  moles  of  external  oxygen  are  present  per  one  mole  of  TNT,  the  detonation 
products  C(s)  and  CO(g)  can  be  oxidized  completely,  according  to  the  following 
reaction: 

3.5C(s)  +  3.5CO(g)  +5.2502(g)  — >  lC02(g)  +  2.373  MJ / moleTNT  (21) 

where  (s)  and  (g)  arc  used  to  denote  the  solid  and  the  gas  phase,  respectively. 

This  reaction  typically  takes  place  in  two  steps:  (i)  An  oxidation  of  C(s)  into  CO(g) 
and  (ii)  a  subsequent  oxidation  of  CO(g)  into  C02(g).  A  comparison  of  Eq.  (20)  and  Eq. 
(21)  shows  that  the  total  molar  heat  of  combustion  of  TNT  (1.034+2.370  = 
3.407MJ/mole  TNT)  is  3.295  times  larger  than  the  molar  heat  of  detonation 
(1.034MJ/mole  TNT).  When  TNT  undergoes  only  a  detonation  reaction  as  defined  by 
Eq.(20),  a  cloud  of  finely  dispersed  soot  is  typically  formed.  In  sharp  contrast,  the 
combustion  of  carbon  and  carbon  monoxide  in  accordance  with  Eq.  (21)  gives  rise  to  a 
fireball  which  expands  with  the  blast  front. 
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It  is  well  established  [22]  that  bulk  water  or  water  mist  can  be  very  effective  in 
hampering  or  even  preventing  the  combustion  reaction  as  given  by  Eq.  (21).  Typically, 
no  fireball  is  formed  during  explosion  while  air  remains  clean  and  water  becomes  sooty. 

While  the  mechanism  by  which  water  hampers/suppresses  the  combustion  reaction  is 
not  well  understood,  two  possible  explanations  are  often  cited: 

(a)  As  demonstrated  earlier,  water  evaporation  can  significantly  decrease  the 
temperature  of  the  detonation-product  gas  mixture,  and,  hence,  lower  the  partial 
pressure  of  CO(g)  reducing  the  thermodynamic  driving  force  for  the  oxidation 
reaction  of  CO(g)  to  C02(g);  and 

(b)  The  aforementioned  reduction  in  temperature  can  lower  the  reaction  rate  constant 
and,  thus,  slow  down  the  kinetics  of  the  combustion  reaction  as  defined  by  Eq.(21). 

A  simple  quantitative  analysis  of  these  two  effects  is  presented  in  the  next  section. 

3.4.2  A  Chemical  Thermodynamics  and  Kinetics  Analyses 

As  explained  earlier  if  one  kilogram  of  water  is  evaporated  for  every  kilogram  of  TNT 
detonated  the  temperature  of  the  resulting  detonation- product  gas-mixture  would  be 
reduced  by  about  one-third  (from  -3500 K  to  -2300K).  A  similar  reduction  in  the  partial 
pressure  of  CO(g)  is  expected  which  would  reduce  the  thermodynamic  driving  force  for 
oxidation  of  CO(g)  into  C02(g)  by  cne-third.  This  finding  suggests  that  the  water- 
induced  reduction  of  thermodynamic  driving  force  for  the  combustion  of  the  detonation 
products  most  likely  plays  a  minor  role  in  the  overall  explosion  mitigation  effects  of 
water. 

In  addition  to  reducing  the  hermodynamic  driving  force  for  combustion,  water- 
induced  reduction  in  temperature  of  the  detonation-product  gas-mixture  also  reduces  the 
rate  of  the  combustion  reaction.  The  forward  rate  constant  for  a  reactbn  such  as  the  one 
shown  in  Eq.  (21)  can  be  defined  by  the  following  Arrhenius  type  relation: 

kf  =A2r|3exp(  — - J-)  (22) 

K1 

where  the  pre- exponential  term  A 2,  the  temperature  exponent  (3  and  the  activation  energy, 
E,  are  characteristic  parameters  for  a  given  chemical  reaction,  R  (=  8.314  J/rnol  K)  is  the 
universal  gas  constant  and  T  the  absolute  temperature. 

The  two  elementary  reactions  associated  with  the  combustion  of  the  detonation- 
product  gas -mixture  are: 


C(s)  +  0.5O2(g)-+CO(g) 

(23) 

CO(g)  +  0.5O2(g)^CO2(g ) 

(24) 

The  temperature  dependencies  of  the  forward  reaction  rate  constants  for  the  oxidation 
reactions  given  by  Eqs.  (23)  and  (24)  are  calculated  using  the  data  from  Ref.[23]  and 
shown  in  Fig.  10.  It  should  be  noted  that  a  logarithmic  scale  is  used  along  the  vertical 
axis  in  Fig. 10.  The  results  displayed  in  Fig. 10  show  that  the  water-induced  reduction  in 
temperature  of  the  detonation-product  gas  mixture  from  -3500K  to~2300K  gives  rise  to 
the  corresponding  reduction  in  the  forward  reaction  rate  constants  of  75%  and  96%  for 
the  reactions  defined  by  the  Eqs. (23)  and  (24),  respectively.  These  findings  indicate  that 


A  COMPUTATIONAL  INVESTIGATION  of  VARIOUS  . . . 

25 

the  water -induced  temperature  decrease  in  the  detonation- product  gas -mixture  can 
substantially  reduce  the  rate  of  the  combustion  reaction. 

It  should  be  noted  that  both  the  oxidation  of  C(s)  given  by  Eq.  (23)  and  the  oxidation 
of  CO(g)  given  by  Eq.(24)  are  associated  with  an  “ ignition  temperature” ,  i.e.,  a 
minimum  temperature  necessary  for  the  combustion  reaction  to  initiate.  The  ignition 
temperatures  for  C(s)  and  CO(g)  are  generally  taken  as  1 173K  and  917K,  respectively 
[24].  An  extension  of  the  thermo-chemical  analysis  presented  in  Section  3.1.2,  and  the 
use  of  the  temperature-dependent  specific  heats  displayed  in  Fig.  3  revealed  that  if  2kgs 
of  water  are  evaporated  for  each  kilogram  of  TNT  detonated,  the  temperature  of  the 
resulting  detonation -product  gas -mixture  would  be  around  900K  giving  rise  to  a 
complete  suppression  of  the  combustion  reaction. 

It  should  also  be  noted  that  the  chemical-kinetics  based  analysis  presented  above  is 
valid  under  the  condition  when  the  combustion  reaction  is  controlled  by  the  rates  of  the 
associated  oxidation  reactions.  When  the  rate  of  the  combustion  reaction  is  controlled  by 
the  mass  transport  of  the  fuel  (C(s)  and  CO(g))  and  oxygen,  the  presence  of  water 
droplets  in  the  combustion -product  gas -mixture  can  affect  the  extent  of  turbulence  and, 
in  turn,  increase  the  rate  of  the  combustion  reaction.  This  phenomenon  was  discussed  in 
Section  3.1.2  and  is  not  considered  in  greater  details  in  this  portion  of  the  work.  As 
stated  in  Section  3.1.2,  the  effect  of  water  droplets  on  the  extent  of  turbulence  is 
considered  to  be  of  a  second  order. 


Fig. 10.  Temperature  variation  of  the  logarithm  of  the  rate  constants  for  the 
forward  reactions  defined  by  Eqs.  (20)  and  (21). 


4.  CONCLUSIONS 


Based  on  the  results  obtained  in  the  present  work,  the  following  main  conclusions  can  be 
drawn: 
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1.  Among  the  four  water-induced  explosive- mitigation  mechanisms  studied  in  the 
present  work,  water  evaporation  which  absorbs  the  major  portion  of  the  thermal 
energy  carried  by  the  detonation -product  appears  to  be  the  dominant  mechanism. 
Water  evaporation  reduces  the  temperature  and  pressure  of  the  detonation-product 
gases  and  can  prevent  their  combustion. 

2.  The  presence  of  water  vapor  in  the  gas  phase  has  also  an  explosion -mitigation 
effect  through  the  potential  reduction  of  the  thermodynamic  driving  force  for  the 
combustion  reaction. 

3.  The  transfer  of  momentum  from  the  detonation-product  gas-mixture  to  water  is  a 
key  phenomenon  controlling  the  break-up  of  water  into  micron-size  droplets. 
Consequently,  this  process  acts  as  a  prerequisite  for  an  effective  explosion- 
mitigation  via  water  evaporation. 

4.  While  the  role  of  the  reduced  shock  speed  caused  by  the  presence  of  water  droplets 
in  the  detonation-product  gas -mixture  in  explosion  mitigation  was  not  analyzed 
quantitatively,  it  appears  to  be  of  a  second  order. 
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Appendix  A:  Rate  of  Evaporation  of  Water  Droplets 

In  this  section,  a  simple  model  is  developed  for  an  assessment  of  the  rate  at  which  water 
droplets  evaporate  when  dispersed  within  hot  detonation -product  gas -mixture.  The 
model  developed  is  based  on  the  following  simplifications  and  assumptions: 

(a)  Water  droplets  are  spherical  in  shape  and  all  have  the  same  radius,  R; 

(b)  Water  droplets  are  distributed  uniformly  within  the  detonation- product  gas -mixture 
so  that  it  can  be  assumed  that  each  single  water  droplet  resides  within  a  cube  with 
an  edge  length,  I.  The  cubes  are  distributed  uniformly  in  space.  Then  each  cube  is 

replaced  with  a  sphere  of  an  equal  volume  with  a  radius,  a  =  (3  / )1/3  /  .  Since 
each  such  spherical  cell  is  identical,  there  is  no  heat  flux  across  its  surfaces. 

(c)  The  volume  fraction  of  the  water  in  the  gas  mixture  can  be  then  defined 

as  f  h2o  =  (R  l  a)3  \  and 

(d)  Due  to  the  intrinsic  spherical  symmetry  of  the  problem  as  posed  above,  a  one¬ 
dimensional  mathematical  model  within  the  spherical  coordinate  system  can  be 
used. 

The  problem  defined  above  can  be  mathematically  expressed  by  the  following  energy 
conservation  equation: 


01  _  A  01 

pgas  p’gas~dt~  ^  7a7+ 


(A.1) 


subject  to  the  initial  conditions: 
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T(r>  R,t  =  0)  =  Tinit'gas  (A.2) 

T(r<R,t  =  0)  =  TinitHiO  (A3) 

and  the  boundary  conditions: 

=  0  (A.4) 

r=a 

'  ^h2o  ~  Ph2o^p,h2oTh,o^h2o  f°r  Th2o  <r^H2o,boii^  and  (A.5) 

r=R 

'  A-H20  —  P H20^-‘H20^ H20^0X  Th20  ~r^H20,boil->  (A.6) 

r=R 

where,  p  is  density,  Cp  constant-pressure  specific  heat,  k  thermal  conductivity,  T 
temperature,  r  spatial  coordinate,  L  latent  heat  of  evaporation,  A  surface  area,  V  volume, 
subscripts  “gas”  and  “H2O”  are  used  to  denote  the  quantities  pertaining  to  gas  and  water 
droplets,  respectively,  while  a  raised  dot  is  used  to  denote  the  first  time  derivative  of  a 
quantity. 

The  boundary  condition  defined  by  Eq.  (A.5)  is  used  while  water  droplets  are  heated 
up  to  their  boiling  point,  TH  0  hojl ,  while  Eq.  (A.6)  is  used  for  the  case  when  water 

droplets  are  at  THOboil  and  evaporation  occurs. 

Eqs.  (A.  1 )  —  (A.6)  are  solved  using  an  explicit  finite -difference  procedure.  A  standard 
mesh  convergence  analysis  is  carried  out  to  ensure  that  the  effect  of  the  mesh  size  has 
been  essentially  eliminated.  The  results  of  this  analysis,  however,  will  not  be  presented 
here  for  brevity. 


Appendix  B:  Non-equilibrium  Liquid  Droplet  Evaporation  Model 


In  a  recent  paper.  Miller  et  al.  [25]  carried  out  a  comprehensive  overview  of  the  existing 
non-equilibrium  two-phase  models  for  evaporation  of  a  liquid  phase  dispersed  in  the 
form  of  single- species  spherical  droplets.  The  droplets  exchange  their  momentum  with 
the  surrounding  carrier  gas  only  via  the  drag  forces,  while  the  thermal  energy  is 
exchanged  only  via  the  convective  heat  transfer.  The  models  had  to  be  modified  before 
they  can  be  used  in  the  present  case  in  order  to  include  the  assumption  that  the  droplets 
are  stationary  relative  to  surrounding  gas.  Under  these  conditions,  evaporation  of  a 


simple  droplet  is  defined  by  the  droplet  temperature,  7]/,  and  the  droplet  mass ,md  , 


evolution  equations  as: 

dL 

dt 


f2Nu 


3Prr; 


f  0!  1 

Tc-Td)+ 

f  L  ' 

lCiv 

in., 


-H, 


(B.l) 


and 
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dmd 

dt 


Sh 


3Scr 


m. 


vTrf  J 


W, 


(B.2) 


where  m  ,  = 


dmd 

dt 


(negative  for  evaporation),  Tc  is  the  local  carrier-gas  temperature,  L 


the  latent  heat  of  evaporation,  0j  —  CpG/CL  is  the  ratio  of  the  constant-pressure  gas 


heat  capacity,  Cp  G  ,  to  that  of  the  liquid  phase  Q  ,  the  gas -phase  Prandtl, 
PrG  —\lCpG/XG,  and  Schmidt,  ScG  =  |dG  /  p GDG  ,  numbers,  (X  gas-phase  viscosity,  A, 

gas -phase  thermal  conductivity,  D  binary  gas-phase  diffusion  coefficient  and  pG  gas- 
phase  density.  The  subscripts  denote  the  droplet  (d),  gas -phase  property  away  from  the 
droplet  surface  (G),  vapor  phase  C),  and  liquid  phase  (L).  In  Eqs. (B.l)  and  (B.2), 
Zd=p dd1 2 3 4 5 6 7 8/(18\vG)  is  the  particle  time  constant  for  Stokes  flow,  where  d  is  the  droplet 
diameter,  and  f2  is  a  heat  transfer  correction  due  to  evaporation,  and  the  Nusselt  (Nu)  and 
Sherwood  (Sh)  numbers  are  empirically  modified  for  convective  corrections  to  heat  and 
mass  transfer,  respectively.  Finally,  H?T  accounts  for  all  additional  terms  used  to 
incorporate  non-uniform  internal  temperature  effects  (i.e.  finite  liquid  thermal 
conductivity),  and  HM  represents  the  specific  driving  potential  for  mass  transfer 
(analogous  to  T(,-Td  for  heat  transfer). 

Eqs.  (B.l)  and  (B.2)  have  been  cast  into  these  specific  forms  in  order  to  highlight  the 
fact  that  various  models  proposed  in  the  literature  differ  predominantly  with  respect  to 
how  the  parameters  f2,  HM  and  H?  arc  calculated.  Five  different  models  (denoted  as  M3- 
M7  in  Ref.  [25])  were  implemented  into  present  work  and  using  the  same  numerical 
procedure  described  in  Appendix  A. 
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